function vmd_trapped(fbase, fnum, OmegaBounce)
%
%   vmd_trapped(fbase, fnum, OmegaBounce)
%
%   Reads file with density distribution.
%   Creates one images.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Domain information and colormap limits are edited in the file.
%   (OmegaBounce = 3.8 at t= 100 in the 2D test run)
%   (OmegaBounce = 5.4 at t= 500 in the 2D test run)
%   (OmegaBounce = 4.1 at t=1000 in the 2D test run)
%   (OmegaBounce = 5.3 at t=2000 in the 2D test run)
%
%   (OmegaBounce = 4.6 at t=1000 in the 3D test run)
%   (OmegaBounce = 4.6 at t=2000 in the 3D test run)
%   (OmegaBounce = 4.7 at t=3000 in the 3D test run)


%---------------------------------------------------------

% system size

jx=48; jy=48; jz=16; jv=256; jdir=1;
%jx=48; jy=1; jz=16; jv=256; jdir=2;
%jx=1; jy=8; jz=32; jv=256; jdir=2;

% physics parameters

vmax=10; Lx=450; Ly=450; 
kz=0.34; nz=1; vPhase=3.53235; phi00=0.01; SpeckleWidth = 99.91139;

%vmax=5;kz=0.34; nz=1; vPhase=3.5; phi00=0.005; SpeckleWidth = 1000;

% colormap bounds

fzvmax = 0.006;


%---------------------------------------------------------

  dv = 2*vmax/(jv-1);
  v  = (0:jv-1)*dv - vmax - vPhase;
  f0 = ( 1/sqrt(2*pi) )*exp( - 0.5*(v + vPhase).**2 );

  dz = 2*pi*nz/(jz*kz);
  z  = (0:jz-1)*dz;

  ixc = min(jx/2+1, jx);
  iyc = min(jy/2+1, jy);


%---------------------------------------------------------
%-- reading density distribution

  fn = num2str(fnum,'%04d');

  fname = [fbase, '.dfu.', fn]

  fid = fopen(fname, 'rb');
  dfu = fread(fid, jv*jx*jy*jz*jdir, 'double');
  fclose(fid);

  dfu = reshape(dfu, jv,jz,jy,jx,jdir);
  dfu = sum(dfu, 5);

  dfv = sum(sum(sum( dfu(:,:,:,:), 4),3),2) / (jx*jy*jz);

  %Ekin=sum(dfv'.*v.^2)*0.5*dv

  fname = [fbase, '.dfv.', fn, '.txt'];
  fid = fopen(fname, 'wt');
  fprintf(fid, '#1.v  2.df(v)\n');
  for i=1:jv 
     fprintf(fid, "%14.6e %14.6e\n", v(i), dfv(i) );
  end
  fclose(fid);

  figure(1, "Name", "averaged df(v)");
  plot(v, dfv, 'o-');
  xlabel("v");
  ylabel("df(v)")

  figure(2, "Name", "averaged v^2 df(v)");
  plot(v, dfv .* v'.^2, 'o-');
  xlabel("v");
  ylabel("v^2 df(v)")


%-- plot fzv --

  fzv =  dfu(:,:,iyc,ixc);

  fzv =  fzv + meshgrid(f0, 1:jz)';

  vlimit = sqrt(2*OmegaBounce /SpeckleWidth/0.37233 /kz);

 
  indv = find( abs(v) < vlimit );

  fzv = fzv(indv, :);

  printf("trapped limits :  [%8.4f, %8.4f]\n",  min(min(fzv(:,:))), max(max(fzv(:,:))) )


  figure(3, "Name", "fzv trapped")
  imagesc(flipud(fzv))
  one_image(fzv, 0, fzvmax, [fbase, '.fvz.', fn, '.png' ])



%---------------------------------------------------------
